Journal of Computational Chemistry

We introduce Bravais Lattice Sampling (BLS), a two-phase method for detecting connected high-density regions in three-dimensional space. BLS places probe sites on a Bravais lattice scaled to the expected nearest-neighbour distance dNN of the target structures, then recovers cluster boundaries by depth-first expansion seeded only from occupied probes, replacing the exhaustive raster scan that conventional connected-component labelling uses to discover seeds. The spacing between probe sites is set from the covering radius of the lattice, which is what allows the method to state in advance the size below which a cluster may escape detection. The second phase, an expansion refinement activated only on probes that return an occupied voxel, verifies every edge, so the components returned are true connected components. BLS versatility allows for selection of different Bravais lattice unit cells to match the target structure; for amorphous, non-crystalline shapes, BLS can default to a simple face-centred cubic unit cell, where the expected minimum cluster size is the only parameter that needs to be set. The current BLS implementation has been developed as a post-processing tool for molecular dynamics trajectories, and was tested for searching water ice clusters of different morphologies. BLS returns component counts and maximum cluster sizes identical to exhaustive-labeller algorithms, with 100% recall; it runs at about 0.94 times the cost of depth-first search, and at 0.84 to 0.90 times the cost of the fastest other labeller in our benchmark set. This algorithm, although implemented by us for molecular dynamics applications, could be of interest in other domain areas where searching for high-density elements in 3D space is relevant.

This work presents DSSNA (Define Spatial Structure of Nucleic Acids) v2026, a free and open-source standalone GROMACS module for the automated analysis of nucleic acid structure based on atomic coordinates. DSSNA reproduces the core functionality of the X3DNA and DSSR approach while extending it to the analysis of molecular dynamics trajectories through integration with the GROMACS software package. The algorithm identifies and classifies structural features, including canonical and noncanonical base pairs, hydrogen-bond networks, stacking interactions, helices, and many other elements of nucleic acid secondary and tertiary structure. Its modular output system enables users to select the structural metrics required for a particular analysis, while the use of neighbor-search algorithms and data-processing capabilities facilitates the analysis of large trajectory datasets. The performance and applicability of DSSNA were demonstrated using a molecular dynamics trajectory of a GUAA tetraloop mutant of the sarcin-ricin domain from Escherichia coli 23S ribosomal RNA (PDB ID: 1MSY). The analysis showed that only a subset of the observed base pairs and stacking interactions remains stable throughout the trajectory. Persistent canonical base pairs formed the structural core responsible for maintaining the secondary-structure profile, whereas most noncanonical base pairs and stacking interactions were transient. The results demonstrate that DSSNA can be used to characterize both static nucleic acid structures and the temporal evolution of structural interactions in molecular dynamics simulations. The work also provides a glossary of terms relevant to the analysis and interpretation of nucleic acid structures.

Nucleic acid secondary structure thermodynamic models have great utility for the analysis and design of nucleic acid structures, devices, and systems, yet their complexity and piecemeal development over decades pose challenges to further innovation. Employing different functional forms for each loop type within a structure, existing nearest neighbor (NN) models necessitate numerous loop-specific dynamic programming recursions in order to calculate physical quantities of interest. Here, we derive and implement a new class of models based on transfer matrices, yielding a single functional form for the free energy of all loop types (as the trace of a matrix product) and a single dynamic programming recursion for calculating physical quantities of interest. This new class of models has additional attractive properties: 1) incorporating a subensemble of states describing proximal interactions within each loop, 2) approximating logarithmic entropic scaling for all closed loops, 3) offering the opportunity to model correlations between the subensembles of neighboring loops by generalizing transfer matrices to tensors, and 4) providing a basis for systematically adjusting the number of model parameters to enhance accuracy as new empirical data become available for parameter regression. We use literature experimental data to regress parameters for representative new transfer matrix models, observing accuracy as good as or better than that of existing NN models, yet using dramatically fewer parameters (e.g., dozens vs hundreds or thousands). For complex ensembles containing an arbitrary number of interacting nucleic acid strands totaling N nucleotides, the new single recursion leads to exact dynamic programs for calculating diverse physical quantities (including the partition function, equilibrium base-pairing probabilities, minimum free energy and suboptimal proxy structures, Boltzmann-sampled secondary structures, and gradient of the partition function) with O(N3) time complexity. The simplicity of the new partition function dynamic program enables overflow-safe partition function calculations using renormalization to achieve scalability for large complexes (e.g., 50,000 nt) as well as formulation using a global matrix multiplication primitive that facilitates GPU acceleration, yielding speedups of an order of magnitude relative to NN partition function algorithms for large complexes. The new class of transfer matrix and tensor models provides an elegant framework for future efforts to advance the accuracy, speed, and scope of computational nucleic acid analysis and design.

Carotenoids are structurally diverse isoprenoid pigments that play central roles in photosynthesis, photoprotection, membrane organization, and cellular signaling. Despite their biological and technological importance, atomistic simulations of carotenoids remain limited by the lack of a transferable force field spanning the chemical diversity of naturally occurring compounds, including glycosylated and acylated derivatives. Here we present CARS (CARotenoidS), a transferable force field for carotenoids that integrates seamlessly with the AMBER family of biomolecular force fields. Parameters were systematically optimized against 22 957 r2SCAN-3c reference energies for 25 representative molecular fragments, yielding an accurate description of polyene conformational energetics, ring rotations, and molecular geometries. Across diverse validation sets, CARS substantially outperforms GAFF2 and provides improved agreement with quantum-chemical reference data for glycosylated and acylated carotenoids. For zeaxanthin, CARS also surpasses OPLS-AA, CGenFF, and previously published carotenoid-specific parameters set in reproducing conformational energetics and structural properties. Two complementary parameter sets are provided: CARS for glycosylated and non-lipidated carotenoids, and CARS+Lipid21 for carotenoids containing saturated or monounsaturated lipid chains. By providing the first unified and transferable parameterization covering the structural diversity of natural carotenoids while remaining fully compatible with established AMBER force fields, CARS removes the need for molecule-specific reparameterization and enables reliable molecular simulations of carotenoids in proteins, membranes, and other complex biological assemblies.

Adaptive sampling accelerates the exploration of conformational space in molecular dynamics (MD) simulations by repeatedly analyzing the accumulated trajectories and seeding a new round of simulations from informative configurations. A growing collection of adaptive sampling policies has been proposed, each built around a particular notion of what makes a configuration informative, yet these methods are scattered across separate and often incompatible implementations, which complicates their systematic comparison and their combined use in meta adaptive sampling schemes. Here, we present AdaptivePy, a compact and extensible Python framework that implements nine seed-selection policies behind a single configuration-driven interface, spanning simple population-based baselines, several established machine-learning and geometry-based methods, and two ensemble or meta sampling policies introduced in this work. We show that the shared implementation reproduces the characteristic selection behavior of each policy on a series of analytic benchmark landscapes. We also introduce a new adaptive sampling scheme that employs TS-DAR, a deep learning framework originally designed to identify transition states, into an acquisition criterion that drives the discovery of an entire multi-basin landscape starting from a single basin. We further demonstrate that the common interface enables meta adaptive sampling policies, which aggregate the rankings of several policies into a single set of seeds. AdaptivePy thereby provides a unified testbed for the adoption, benchmarking, and continued development of adaptive sampling methods for biomolecular MD simulations.

Aptamers are structured nucleic acid ligands capable of high affinity, high specificity molecular recognition generated using variations of the SELEX (Systematic Evolution of Ligands by Exponential Enrichment) process. However, SELEX often produces sequences that enrich yet may lack binding efficacy. We propose a measure called the Ruggedness Composite Index (RCI) along with a method for computing it, that can be used to distinguish binding-competent ('active') aptamers from weak or non-binding ('inactive') aptamers. Given a set of aptamers, RCI incorporates information on their fragmentation (landscape partitioning), basin entropy (metastable state distribution), cumulative density irregularity (non-uniform occupancy), and structural energy correlation length (structure-energy coupling scale). We test whether secondary-structure folding energy landscape topology distinguishes active from inactive aptamers using a multiscale level set framework across six datasets. Active aptamers show lower RCI values and occupy smoother, funnel-like conformational spaces, while inactive aptamers show higher RCI values, reflecting fragmented, high-entropy landscapes. By contrast, classical thermodynamic features, such as minimum free energy, show limited discrimination between active and inactive aptamers. In all datasets, sequences that exhibit enrichment which is not monotonic but lack specificity exhibit elevated ruggedness, indicating landscape topology can predict non-specific enrichment. These results indicate that folding landscape organization can be used as a predictor of aptamer activity and establish RCI as a simple, mechanistically interpretable measure for improving candidate prioritization, especially in therapeutic aptamer discovery.

Oncogenic mutations are known to occur at non-random rates and specific hotspots, yet the structural and thermodynamic factors underlying these hotspots remain poorly understood. This study presents molecular dynamics simulations of 22 cancer hotspots across 10 oncogene families, simulated as histone-DNA complexes. Across 18 of 22 mutations, thermodynamic effects were predominantly localized to within 6 [A] of the DNA-histone interface (mean capture 104%), with van der Waals interactions driving the effect at the thermodynamic extremes, indicating that oncogenic mutations alter nucleosome binding through precise local contact changes rather than global structural rearrangements. Further, analysis revealed a gene-family correlated pattern in thermodynamic stability. RAS family mutations showed a consistent trend toward nucleosome stabilization (mean {Delta}{Delta}G = -4.82 kcal mol-1), while kinase domain mutations trended toward destabilization (mean {Delta}{Delta}G = +40.12 kcal mol-1), a difference reaching statistical significance in this exploratory analysis (Mann-Whitney U, p = 0.008). This interface-specific mechanism, combined with the gene-family-correlated thermodynamic pattern, provides a biophysical framework for understanding nucleosome-level contributions to cancer hotspot mutation biology.

We used simulations of short single stranded DNA (ssDNA) homopolymers, based on the sequence dependent Three Interaction Site (TIS) model, to calculate the persistence length (lp) in polyvalent cations. The TIS model accounts for stacking interactions and electrostatic interactions are treated using the Coulomb potential. We find that lp for dT30 (T is thymine) and dA30 (A is adenine) is quantitatively fit using [Formula] ([Formula] is the bare persistence length,{lambda} is a dimensionless constant, and{kappa} -1 is the inverse Debye length) in the divalent cations Mg2+ and Ca2+. The dependence of lp on{kappa} is surprising because it was derived for long flexible polyelectrolytes in which the charges interact via the Debye-Huckel potential. The [Formula] values are 0.4 nm and 1.1 nm for polyT and polyA, respectively. Strikingly, lp is almost independent of the tetravalent spermine concentration. There is no clear theoretical explanation although simulations suggest that the number of spermine molecules that bind to the ssDNA saturates at a small value. A qualitative picture, based on the restrictions of access to the phosphate groups due to volume exclusion of the anisotropic structure of Spm4+, rationalizes the simulation results. The predicted dependence of lp in spermine awaits experimental test.

Interleukin-13 (IL-13) is an immunomodulatory cell signaling cytokine that has been implicated in neurodegenerative disease and chronic inflammation. IL-13 binds to its low and high affinity receptors, IL-13 receptor 1 (IL-13R1) and IL-13 receptor 2 (IL-13R2), respectively, with residence times that vary accordingly. As the binding kinetics of the cytokine-receptor complexes influence cellular responses, we employed the molecular dynamics (MD) simulation-based{tau} -random acceleration molecular dynamics method ({tau}RAMD) to compute relative residence times for wild-type (WT) IL-13 and 19 IL-13 mutants to the two receptors. Comparison with experimental kinetic data shows that the{tau} RAMD computations capture the trends in residence times. Analysis of simulated dissociation trajectories of the cytokine-receptor complexes reveals two distinct dissociation pathways of IL-13 from each of the receptors. This study thus pinpoints key determinants of the interaction of IL-13 with its receptors which could be targeted for therapeutic design. Statement of SignificanceCytokines are regulatory proteins that bind to cell surface receptors and thereby send signals to the cellular interior. Interleukin-13 (IL-13) is a cytokine that has a low and a high affinity receptor. It has important physiological roles, and its deregulation is involved in diseases such as atopic dermatitis and asthma. We employed a molecular dynamics simulation-based method to compute the effects of changes in the sequence of IL-13 on the lifetimes of complexes of IL-13 and its receptors. Comparison with experiments supports the validity of the computational approach and analysis of the simulations reveals two distinct ways in which IL-13 dissociates from each receptor. These results thus provide a map for targeting IL-13 - receptor interactions for the design of therapeutics.

Atomistic simulations of intrinsically disordered proteins (IDPs) are notoriously sensitive to force field inaccuracies, either regarding the protein or the water model, yielding inaccurate observables such as compactness, secondary structure propensities, or kinetics. The currently most widely used IDP force fields are Amber99sb-disp (A99disp) and Charmm36m (C36m). A99disp includes a new water model and thus optimized both, the protein and the water interactions. In contrast, C36m used the Tip3p water model and optimized only protein interactions. In many cases, C36m+Tip3p underestimates radii of gyration compared to FRET or SAXS experiments. Such overly compact structural ensembles are believed to arise from an imbalance between protein-protein, protein-water, and water-water interactions, which might be due to Tip3p inaccuracies. Here, we aim at re-balancing these interactions by combining C36m with the Optimal Point Charge (OPC) water model. Recently, this C36m+OPC combination showed improved accuracy for the disordered domain of the measles virus nucleoprotein. Here we present a systematic assessment, comparing C36m+OPC, C36m+Tip3p, C36m+Tip4p, C22*, A03ws, A99sb-ws, and A99disp for five IDPs, as well as a subset of those for five globular proteins, a set of disordered AGQ-repeat peptides, and the fast folding miniprotein CLN025. We compared extensive MD simulations (> 8.5 ms) with SAXS, NMR, circular dichroism, photo-induced electron transfer (PET), T-jump infrared spectroscopy, and X-ray crystallography measurements. We found that combining C36m with OPC improved accuracy over C36m+Tip3p for IDP ensembles without compromising its accuracy for globular proteins. While also the kinetics of the AGQ-peptides were more accurate for C36m+OPC, those of CLN025 folding were less accurate. Overall, C36m+OPC showed similar accuracy as A99sb-ws and A99disp, the latter is currently considered among the most accurate protein force fields.

DNA often functions in heterogeneous molecular systems containing proteins, lipids, polymers, and other materials. All-atom molecular dynamics simulations can be used to study DNA in these multicomponent systems, but computational cost limits the accessible system sizes and time scales. Coarse-grained models extend these scales, but existing DNA models are generally not designed for interactions with a broad range of other molecular species. To fill this gap, we develop a coarse-grained model of DNA designed for use with the Martini 3 force field. The model was parameterized through an iterative Bayesian optimization workflow, which used a scaled Wasserstein metric to compare distributions of local geometrical features and global structure from coarse-grained simulations against all-atom reference simulations. The optimized model captures key structural and mechanical properties of single- and double-stranded DNA across varying strand lengths and ionic conditions, while retaining compatibility with the broader Martini 3 ecosystem. This compatibility enables DNA to be integrated with a broad range of molecular systems, as we illustrate through simulations of double-stranded DNA bound to a transcription factor, cholesterol-tagged DNA duplex interacting with a lipid bilayer, a crossover-containing DNA nanostructure, and single-stranded DNA adsorbing onto graphene. Together, these results establish a transferable coarse-grained model of DNA for simulations of heterogeneous biomolecular and engineered systems.

The process of drug unbinding is of immense importance in the field of biophysics and therapeutics. The behavior of these systems is greatly influenced by their thermodynamic and kinetic properties. Therefore, it is crucial to accurately estimate the ligand binding free energies and rate of ligand dissociation, yet these processes are often governed by rare event transitions that lie beyond the reach of standard brute-force molecular dynamics simulations. While enhanced sampling simulations offer a solution, their efficacy is strictly contingent upon the selection of appropriate collective variables (CVs) which is non-trivial for complex systems like protein-ligand complexes. In this study, we present a method to derive optimized CV from transition state region (TS) via an interpretable machine learning (ML) model, Elastic Net. By employing some physically intuitive order parameters, the derived optimized CV from the TS-region greatly accelerate ligand binding-unbinding transitions and achieves rapid free energy surface (FES) convergence across diverse systems including buried and solvent exposed active sites such as Trpsin-benzamidine complex, host-guest systems and sodium epoxidase etc. Intriguingly significant contribution of the ligand hydration is found in the optimized CV which depicts crucial role of solvent in driving ligand binding-unbinding transitions. The estimated binding free energies for different protein-ligand complexes match quite well with experiments, while maintaining a low computational cost. The derived optimized CV is also used to calculate the ligand residence times across different systems and calculated residence times are within the experimental range for all systems, again with very little computational costs. Moreover, we show that the optimized CV constructed from TS region via an interpretable ML model is transferable across diverse systems, offering a robust and scalable framework for drug discovery and investigation of complex biomolecular recognition.

The comparison of protein structures is necessary to determine molecular functions and evolutionary relationships and can also facilitate drug design and the prediction of separation options. The number of available protein structures is rapidly increasing, driven by experimental determination but also prediction methods such as AlphaFold. In this context, traditional structure comparison based on atomic superposition becomes computationally inefficient. To enable scalable analysis and surface shape comparison, alternative approaches represent protein structures using compact, fixed-length vectors known as descriptors. Current vectors often contain hundreds of entries, such as three-dimensional Zernike descriptors composed of 121 entries, or combinations of molecular and geometric properties. Here, we propose to use the Laplace-Beltrami spectrum, a mathematical representation of object surfaces derived from the eigenvalues of the Laplace-Beltrami operator, as an efficient option to capture protein structural properties and facilitate structural classification. Specifically, protein surfaces are encoded using only the first 10 non-zero eigenvalues of this spectrum. We used the SHREC 2025 dataset as a benchmark and achieved 85.8% accuracy on the original 97 protein structure classes and 97.7% accuracy on the homology-grouped 45 classes. Our method is fast, requiring ~75 min computation time for the 11,555 protein surface meshes of the SHREC 2025 dataset on a conventional AMD Ryzen 7 5700X 8-Core processor with 32 GB RAM. The approach can easily be applied to large protein structure databases because it is readily parallelizable and allows the incorporation of other descriptors, including surface features such as charge. This rapid screening tool can therefore be used to identify proteins with related surface shapes independent of sequence homology.

Alzheimer's disease is characterized by the accumulation and aggregation of amyloid-{beta}(A{beta}), but the molecular mechanisms linking environmental and infectious factors to A$\beta$ conformational changes remain incompletely understood. Herpes simplex virus type 1 (HSV-1) has been proposed as a potential contributor to AD pathology, and interactions between the viral glycoprotein B (gB) and A$\beta$ may influence the conformational behaviour of the peptide. Molecular dynamics (MD) simulations provide atomic-scale information on such interactions, but conventional structural descriptors may not fully capture changes in the organization of residue interaction networks. Here, we introduce a graph-geometric framework based on Forman-Ricci curvature to characterize the evolution of residue interaction networks during MD simulations. Each simulation frame is represented as a residue interaction graph based on C--C contacts, and residue-wise curvature profiles are analysed across time. We apply the framework to A{beta}1-42 in isolation and in complex with HSV-1 gB. Conventional MD analyses indicate stable association of the simulated complex, favourable interaction energetics, and conformational changes in A{beta}, including a transition from -helical structure toward {beta}-turn-rich conformations over the simulated timescale. Forman-Ricci curvature reveals pronounced and spatially localized remodelling of the A{beta} residue interaction network in the complex, with the strongest changes concentrated in the C-terminal region. These regions also exhibit reduced temporal curvature fluctuations and progressively distinct geometric behaviour throughout the simulation. Hierarchical clustering further identifies cooperative groups of residues with coordinated curvature dynamics, including a prominent C-terminal domain. Together, these results demonstrate that Forman-Ricci curvature provides a complementary description of biomolecular dynamics by capturing changes in the geometric organization of residue interaction networks that are not directly represented by conventional structural descriptors. The framework provides a general computational approach for studying network-level structural remodelling in protein molecular dynamics and offers a quantitative perspective on the conformational consequences of HSV-1 gB--A{beta} association.

Neutron scattering is a powerful technique to investigate atomic structures and molecular dynamics of proteins at the nano-scale. When it comes to dynamics, incoherent and coherent scattering respectively provide information on the single and collective dynamics of nuclei. In proteins, hydrogen has the highest incoherent cross-section, and it is common practice to overlook the contribution of coherent terms stemming from all nuclei. However, the fast collective dynamics of heavier nuclei could also be studied if coherent scattering and incoherent scattering were experimentally separated. The recent advent of polarized neutron spectroscopy with sufficient flux and energy resolution has made it possible, and opens new perspectives to investigate the relative importance of coherent scattering and the information it provides on biological samples. The present study reports on the use of polarized quasi-elastic neutron scattering (QENS) and the application of a minimalistic model adapted to both individual and collective dynamics. Using a perdeuterated green fluorescent protein as a model globular protein, the study provides an interpretation of the dynamical parameters obtained with QENS, and a comparative study of the Elastic Coherent and Incoherent Scattering Factor. Based on both experiments and calculations, we discuss the relative importance of distinct and self components of coherent scattering, which is often wrongly assumed to be representative of collective dynamics only. The results highlight the current impediments rendering complicated a straightforward analysis of fast collective dynamics in hydrated protein samples.

Umbrella sampling is widely used to estimate biomolecular free energy landscapes and relative folding stabilities. Although equilibration is a critical component of umbrella sampling workflows, the impact of equilibration length on thermodynamic predictions remains poorly understood. Here, we investigate the effect of equilibration length on relative folding free energy predictions for four RNA hairpins with loop sequences GUGAAA, CUGGGA, GUAAUA, and UUAAUU with helical stems of three base pairs. Umbrella sampling simulations were performed using an end-to-end distance reaction coordinate spanning 15-45 [A], where equilibrium simulations (windows) were spaced at roughly 1 [A] intervals. In these calculations, the hairpin stem-loops were allowed to equilibrate in an end-to-end distance window and then the coordinates were transferred to the next larger end-to-end distance window to equilibrate. Two equilibration lengths, 2 ns and 100 ns per window, were followed by 600 ns of production sampling. Potential of mean force (PMF) profiles were reconstructed using the Weighted Histogram Analysis Method (WHAM) and used to calculate pairwise free energy differences with thermodynamic cycles. Increasing the equilibration length produced substantial, sequence-dependent changes in the reconstructed free energy landscapes. The 100 ns protocol generated markedly flatter PMFs for GUGAAA and UUAAUU and pronounced reshaping of the free energy landscape for GUAAUA. These changes were accompanied by reductions in hydrogen-bonding and stacking interactions, particularly within the intermediate regions of the reaction coordinate. The resulting thermodynamic predictions, as free energy change differences, were therefore highly sensitive to equilibration length. Across nearly all hairpin pairs, the 100 ns equilibration yielded substantially larger magnitude free energy change difference values than the corresponding 2 ns equilibration, with differences that greatly exceeded replica-to-replica variability. Comparison with optical melting measurements and nearest-neighbor thermodynamic predictions revealed that 2 ns equilibration times more closely agreed with experimental values than those obtained using 100 ns equilibration. These findings demonstrate that longer equilibration can systematically alter the structural ensembles sampled during umbrella sampling and amplify predicted stability differences without improving agreement with experiment. More widely, our results highlighted equilibration length as a critical and nontrivial parameter in RNA free energy calculations and demonstrate that increased equilibration does not necessarily lead to more accurate thermodynamic predictions.

Biomolecular condensates formed by intrinsically disordered proteins require molecular models that accurately describe proteins in both dilute solution and condensed phases. Explicit-solvent coarse-grained models offer an attractive balance between chemical resolution and computational efficiency. Yet, it remains unclear whether improving dilute-state properties is sufficient to obtain an accurate description of condensates. Here, we address this question by introducing minimal modifications to the Martini 3 force field that combine recent advances in bonded interactions with refined protein-water interactions and strengthened glycine self-interactions, while preserving the underlying chemical transferability of the model. The resulting model substantially improves the description of single-chain conformations across a diverse benchmark of disordered proteins. We then investigate phase separation of the well-characterized low-complexity domain of heterogeneous nuclear ribonucleoprotein A1 and its sequence variants. The model reproduces several key physicochemical properties of biomolecular condensates, including chain expansion in the dense phase, sequence-dependent intermolecular contacts, protein diffusion and its relation to single-chain dimensions, and hydration, while revealing quantitative limitations in condensate density, phase equilibria, and ion partitioning. Our results show that improving dilute-state behaviour translates into a better description of condensed-phase properties, including condensate density, but is not sufficient to quantitatively reproduce the equilibrium between the dilute and dense phases.

Background: RNA molecules form thermodynamic ensembles, but interpretation often requires a single representative structure. Existing base-pair centroid estimators assess agreement at the level of individual base pairs and do not directly target nesting depth along the sequence. Methods: We introduce Mountain Centroid, which minimizes expected squared mountain-profile distance, and derive dynamic programming algorithms with and without RNA pairing constraints. We also combine the Mountain Centroid objective with the base-pair centroid gain. Results: Across 21,254 RNAStrAlign sequences, Mountain Centroid had lower median normalized mean squared mountain distance (NMSMD) than minimum-free-energy (MFE) and base-pair centroid ({gamma} = 1) structures, whereas its median base-pair F1 was lower. Imposing RNA pairing constraints improved base-pair F1 for 59.35% of sequences and reduced it for 3.58%. At an illustrative weight, the combined objective had median base-pair F1 similar to MFE while retaining lower median NMSMD than MFE and all tested {gamma}-centroid settings. Conclusions: Mountain Centroid represents an RNA structural ensemble with a single secondary structure that reflects how nesting depth varies across nucleotide positions. Combining mountain-profile and individual-base-pair criteria allows their relative contributions to be varied.

Many functional RNA motifs adopt structures that deviate from the canonical A-form helix and are emerging targets for RNA-directed therapeutics. The microtubule-associated protein tau (MAPT) A-bulge motif (5'-GCAGU/5'-ACGU) is one such motif. Because its structure is stabilized by a delicate balance of local interactions, its accurate modeling remains a major challenge for molecular dynamics (MD) simulations. The experimentally determined nuclear magnetic resonance (NMR) structure of the MAPT A-bulge motif provides a stringent test of whether RNA force fields can accurately reproduce the experimentally observed conformation. Most current AMBER-family RNA force-field models have incorrectly favored a non-native base-triple state of the MAPT A-bulge motif over the experimentally observed stacked state. Structural comparison of the stacked and base-triple conformations revealed that overly favorable NH-N hydrogen bonds between the bulged adenosine and an adjacent Watson-Crick base pair were the primary source of this imbalance. We developed gHBfix-18Ab, an 18-component hydrogen-bond correction that distinguishes NH and NH2; donors. gHBfix-18Ab was combined with the previously developed OL3CP and NBfix0BPh corrections to generate the composite model gHBfix-18Ab*. This model restored the experimentally observed stacked state as the global minimum in the calculated free-energy profile and improved agreement with NMR-derived distance data for the A-bulge region. Importantly, gHBfix-18Ab* did not produce marked structural destabilization of the cUUCGg tetraloop, a widely used benchmark for RNA force-field validation, suggesting that the refinement preserves the stability of the unrelated RNA motif. These results demonstrate that targeted refinement of hydrogen-bond interactions provides a practical strategy for systematic improvement of RNA force fields toward more accurate modeling of noncanonical RNA motifs.

Circuit Topology (CT) provides a fundamental framework for analysing folded polymer chains, with applications in functional annotation, protein engineering and drug development. We present a protein CT analysis plugin for PyMOL v3.1.6.1 with a graphical user interface (GUI), automatic installation, and novel features developed through integration with PyMOL's application programming interface (API). The plugin integrates various previously developed CT methodologies for studying structured proteins and their complexes as well as the dynamics of disordered proteins. Analysis of a representative protein and a molecular dynamics trajectory demonstrates the plugin's three analysis modes and their outputs. The plugin reproduces the reference ProteinCT implementation exactly on the structures tested, and is distributed with a versioned release, a pinned environment and a one-command reproduction of every result reported here.